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ABSTRACT 



A high resolution acoustic navigation system for ocean use 
is being developed at the Woods Hole Oceanographic Institution. 
The system can yield navigation fixes with respect to a bottom 
moored reference net with accuracies (on a fix to fix basis) of 
a few centimeters. In order to use the system to best advan- 
tage a survey is required to determine precisely the relative 
positions of the net elements. Each element combines a pulse 
transponder with a continuous wave (CW) beacon. Accumulated 
phase (Doppler shift of the CW beacon) between survey points is 
measured as well as acoustic travel times between survey points 
and transponders. Non-linear regression techniques are employed 
to develop a maximum likelihood estimator for net element 
positions based on these phase and travel time measurements. 

An approximate error covariance matrix is generated and an 
optimum choice of survey points is indicated. The combined 
system, using these selected locations for performing the survey, 
can yield reference mooring coordinates with error of ±1 meter. 
Improved precision appears to be limited by inaccuracies in the 
pulse and Doppler measuring system. 



Thesis Supervisor: Robert C. Spindel 

Title: Assistant Scientist at Woods Hole Oceanographic 

Institution. 

Thesis Supervisor: Arthur B. Baggeroer 

Title: Associate Professor of Ocean and Electrical Engineering 

at Massachusetts Institute of Technology. 



- 2 - 



ACKNOWLEDGEMENTS 



I would like to express my gratitude to all those who have 
provided assistance and support in this endeavor. In particular 
I would like to acknowledge the following: the U.S. Navy for 
financial support of my graduate education; my advisors, Drs. 

R. C. Spindel and A. B. Baggeroer for their advice and encourage- 
ment; Drs. R. P. Porter and W. K. Smith for their hours of 
discussion, criticism and comments; Gatha Mosier for her valuable 
secretarial assistance; and finally, my wife and daughter for 
their patience and understanding. 



- 3 - 





I 




TABLE OF CONTENTS 

PAGE 

ABSTRACT ' 2 

ACKNOWLEDGEMENTS 3 

LIST OF FIGURES 6 

I. INTRODUCTION 7 

1. Description of System ^ 

2. Previous Survey Techniques H 

3. Pulse-Doppler Survey Technique ^2 

II. PARAMETER ESTIMATION 15 

1. The Maximum Likelihood Estimate 15 

2. Properties of the Maximum Likelihood Estimator* •• • 17 

3. The Maximum Likelihood Estimate in the Case of 

Additive Gaussian Noise 1® 

4. Application of Parameter Estimation Techniques 

to the Survey Problem 22 

5. Gauss-Newton Method of Non-linear Least Squares 

Estimation 26 

6. Modifications of the Basic Gauss-Newton Procedure . 28 

7. Iterative Technique 29 

8. Error Estimates 30 



- 4 - 




1 

\C-. 

I 



i 



1 



I 



PAGE 



III. MINIMIZATION OF ERRORS 33 

1. Measurement Errors (Pulse System) 33 

2. Measurement Errors (Doppler System) 40 

3. Selection of Survey Coordinates 45 

4. Computer Simulation 47 

IV. DISCUSSION 65 

V. APPENDICES 68 

A. Description of Program ARCOR 68 

B. Detection Probabilities of W.H.O.I. Cycle Counter. 85 

C. Sample Calculation of Signal-to-Noise Ratio 93 

BIBLIOGRAPHY 95 



- 5 - 




I 



I 

I 

I 

1 



LIST OF FIGURES 

PAGE 

Figure 1. Acoustic Navigation Beacon Net 8 

Figure 2. Survey Coordinate System 23 

Figure 3. Effective Sound Speed vs. Geometric 

Angle 39 

Figure 4. W.H.O.I. Doppler Receiver (Block Diagram).. 41 

Figures 5A-5E. Arrangement of Survey Locations 48 

Figure 6. Program PDSRV (Block Diagram) 55 

Figure 7. Variables Used in Program ARCOR 69 

Figures 8A-8C. Probability Densities of Signal and Noise 

at W.H.O.I. Cycle Counter 88 



- 6 - 



I 






I INTRODUCTION 



There are many methods currently used to position ships, 
submersibles, buoys or submerged instruments in the deep ocean. 
Broadly speaking they can be divided into two types. The first 
method positions via electromagnetic transmission from a shore 
station or satellite whose position or orbit is precisely known 
Such systems operate at long ranges but are limited in accuracy 
to about 100-200 meters at best. Furthermore, these systems 
require above surface antennas and are therefore of little use 
in positioning entirely submerged vehicles. The second method 
positions via acoustic transmission from a set of transponders 
or beacons, usually moored to the ocean bottom, whose relative 
positions are precisely known. When operated in the pulsed or 
transponder mode, position inaccuracies of an acoustic system 
can be as small as 10-20 meters in 5 km water depths. When 
operated in a continuous, or Doppler mode, the accuracy can 
approach 3-4 centimeters. The range of acoustic systems is 
usually severely limited by a number of factors, the most im- 
portant of which are acoustic refraction and attenuation of 
sound in sea water due to spreading and absorption in the ocean 
Current acoustic systems have ranges of about 10 km. 

A navigation system utilizing acoustic beacons (Fig. 1) 

I 

v/hich operate in both a pulse and continuous wave mode is being 
developed. This system capitalizes on the attributes of both 
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Figure 1. Acoustic Navigation Beacon Net 
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the pulse and Doppler modes and will be capable of positioning 
a platform with respect to an array of reference beacons to an 
accuracy of 1-2 meters and of repositioning a platform within 
10 centimeters of a previous position. To make full use of the 
system capabilities a survey is required to accurately deter- 
mine the relative positions of the ref-erence beacons. The pur- 
pose of this investigation is to develop an accurate survey 
technique. Algorithms for estimating the beacon positions and 
the errors in those positions are derived. Performance is 
predicted using computer models. The estimation method used is 
a straightforward extension of the pulse system technique des- 
cribed by Hunt ^ a^. (1974) . 

1) Description of System 

Pulse (transponder) navigation systems measure the travel 
time of an acoustic pulse to estimate the slant range to a 
transponder. Shipboard processing equipment interrogates the 
bottom moored transponders and measures the time of arrival of 
the reply pulse from each transponder. Sound velocity correc- 
tions and corrections for fixed system delay times are applied^ 
and travel times are converted into estimates of slant range to 
the transponders. The accuracy of the measurement is depen- 
dent upon the pulse bandwidth and the signal-to-noise ratio in 
the received pulse. Systems of this type are widely in use in 
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oceanographic research, submerged rescue operations and commer- 
cial enterprises in the ocean (Boegeman et aj^. , 1972; Cestone 
and St. George, 1972; Van Calcar, 1969). 

A continuous wave (CW) beacon system has been developed 
at the Woods Hole Oceanographic Institution to accurately 
track the position of ship-mounted and 'submerged hydrophones 
(Porter ^t al^. , 1973). The system measures the frequency 
change of a continuous tone due to hydrophone motion. Hydro- 
phone velocity and displacement are derived from the beacon 
signal Doppler shift which is given by 

where fj^ is the Doppler shifted beacon frequency due to a 
hydrophone velocity v, f„ is the beacon frequency and c is the 
speed of sound in water. Doppler shift is proportional to 
beacon frequency. Hence the accuracy of the measurement is 
dependent on the beacon frequency and the signal-to- noise 
ratio in the narrow band that encompasses the maximum Doppler 
shift. Because the beacon frequency is much greater than the 
pulse bandwidth of transponder systems the Doppler system can 
measure range changes much more accurately than the pulse 
system assuming that both systems have the same signal- to - 
noise ratio. In the present implementation of the system^ 
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accumulated phase change 0 {t) = 2 Tr(Af)t is measured by digital 
techniques and hydrophone displacements .of a quarter wavelength 
can be resolved. A wavelength for 13 kHz (a typical beacon 
frequency) is ^ = 11.5 centimeters, so that resolution is 
about A /4 ~ 3 cm. Phase changes due to inhomogeneities and 
fluctuations in the water column are negligible over the range 
of the system. We discuss this in Section III. 2. 

2) Previous Survey Techniques 

Several investigations have been conducted to develop an 
optimum technique for determining the coordinates of ocean 
bottom acoustic transponders. No satisfactory survey has been 
developed previously for an array of bottom-moored CW beacons. 

Surveys of acoustic transponders generally fall into iwo . 
categories: baseline crossing and iterative techniques. The 
baseline crossing technique, described by Haehnle (1967) and 
Hart (1967), is often referred to as the conventional or class- 
ical transponder survey method. The method requires that the 
depth of the transponders be previously detemined from inde- 
pendent measurements. The baseline length, i.e., the horizon- 
tal distance between two transponders^ is determined by steaming 
across a baseline and measuring the slant range to the two 
transponders during the traverse. This requires accurate ship 
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positioning and large amounts of ship time to yield accurate 
results. More recent iterative techniques are based on or are 
similar to that described by Vanderkulk (1961). Vanderkulk 
showed that for three bottom transponders and six co-planar 
survey positions whose depths are known, the positions and 
depths of the transponders can be found by solving six linear 
equations. A unique solution exists when the survey positions 
do not lie on a conic section. When more than six survey posi- 
tions are used the additional equations result in an over- 
determined situation for which minimum least squared error 
techniques have been successfully applied, (Lowenstein, 1965; 
McKeown, 1974; Mourad et , 1972; Heckman and Abbott, 1973; 
Hunt et , 1974) . 

3) Pulse-Doppler Survey Technique 

The pulse-Doppler survey consists of making measurements 
of travel times of acoustic pulses and accumulated phase 
changes of CW signals at various survey points. These are con- 
verted to slant ranges and slant range differences, respect- 
ively. The pulse-Doppler survey can yield more accurate results 
than the pulse survey alone due to the higher resolution of the 
Doppler system measurement. 

The slant range (SR) between a transponder and a survey 
point is equal to the travel time (TT) multiplied by an 
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effective sound speed, Ve: 



SR = Ve • TT 

The effective sound speed is a function of the sound velocity- 
profile, the depths of the transponder^ and the survey point 
and varies with the amount of acoustic refraction (Vass, 1966) . 
For m survey points there are m measurements of slant range 
for each transponder. 

The change in slant range, DSR, between a beacon and two 
survey points is proportional to the change in accumulated 
phase, 0 (.^ 2 ) - 0^(t^^) and the wavelength of the signal, \ = c/f^: 

DSR = [ eT(t2) - !2^(tj^)] X /2 TT . 

The received signal, 

s(t) = A(t) exp { -j2 TT f^-t } , 
where fj^ = fg +fgv/c is the Doppler shifted frequency^ can 
be written as 

s(t) = A(t) exp { -j [2 TT f^t + gr(t) ] } 
where the phase variation due to a finite receiver velocity is 
0 {t) = 2 IT fgt v/c (Porter ^ £1^., 1973). In general v is a 
function of time, v = v(t). However, we assume that the 
sampling interval is short enough to consider v constant. in 
practice, an interval of .1 sec is used so the assumption is 
valid. Measurement of accumulated phase change provides 
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velocity and displacement information along the beacon-survey 
point vector, i.e., the displacement is. the change in slant 
range. For m survey points there are m - 1 measurements of 
slant range differences for each beacon. 

The measurements of travel times and phase changes may 
be made at the same or different survey points. Only the case 
where the measurements are made at the same points will be 
discussed in detail since this technique is more efficient in 
practice. The method of solution for the transponder coordin- 
ates is the same in either case but the associated errors will 
be different. 

Non-linear regression methods as outlined by Draper and 
Smith (1966) are used to analyze this survey problem. Under 
the assumptions given in the following section the maximum 
likelihood estimator for the beacon positions is found. An 
approximate covariance matrix for the estimates and a pro- 
cedure for optimizing the choice of survey points is developed. 
The sources of error and their effect on the estimates of the 
beacon coordinates are discussed. 
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II PARAMETER ESTIMATION 



The problem of parameter estimation can be briefly des- 
cribed as follows: given a system from which information can 
be received and whose state is characterized by a set of p 
parameters, make some estimate of the state of the system from 
the information received. The information, which will gener- 
ally be perturbed by noise, may be measurements, a message, or 
some other observations on the system. 

The following definitions are used: 

= [9]^/ 02 ' = true parameter values 
= [0j^, 02 »'*'/ 0 p] = estimated parameter values 

F' = F ' (0) = [fj^, f 2 , . . . fj^ ] = noiseless message 

Y' = [y^^, y 2 » - ' • » yj^ ] = message (observations perturbed 

by noise) 

where 9 _, 9 ^, and Y denote column vectors and ( ' ) denotes 
transpose (Manasse, 1960) . 

When a noiseless message F(0)is received, the information 
is sufficient to determine exactly if the number of measure- 
ments n is not less than the number of parameters p. When 
n = p the problem is called the MINIMUM DATA case. When n > p 
it is called the REDUNDANT DATA case. 
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1) The Maximum Likelihood Estimate 



When the observations Y are received, some criteria must 
be used to determine an optimum estimate ^ since Y differs 
from ^ due to the presence of noise. One often used decision 
rule is based on the maximum likelihood principle. This 

A 

principle prescribes that the observer choose the £ which 
renders the observations Y most likely; that is, it maximizes 
the conditional probability density P ( Y | ^) . When expressed 
as a function of Y and ^ this is called the likelihood function 
and is sometimes written L(Y;9). When the ^ priori probability 
density P (^) can be considered to be a constant over the 
region of interest then the maximum likelihood estimator also 
maximizes P(^|y). This can be shown by the use of Bayes 
Theorem to express the ^ posteriori probability density 
P I Y) in terms of P ( Y | 0^) and P (9 ) : 

P(y|9) P(9) 

P(9|Y) = = C3^P(Y|9) (1) 

P(Y) 

where P ( Y) = / P(y|£) P(9) d9 

(Note that P(Y) is not a function of ^ and can be treated as a 
constant of proportionality) . 

The estimate ^ which maximizes the likelihood function 
L(Y;^) must satisfy the equations 
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Small variance is a desired property. An asymptotic measure 
of this property is efficiency. An estimator is said to be 
efficient if the asymptotic variance of the estimate is no 
larger than the asymptotic variance obtained using any other 
estimator. The maximum likelihood estimator, when it is con- 
tinuous and the first and second derivatives exist and are 
absolutely integrable, is efficient. It can be shown that 
maximum likelihood estimates are asymptotically normally dis- 
tributed with mean ^ and variance 1/R (^) where 



- f 

In the case of additive Gaussian noise the form of the 

I likelihood function can be determined and expressions for the 

I 

, variance of the estimates can be developed. 

' 3) The Maximum Likelihood Estimate in the Case of Additive 

Gaussian Noise 

The message is assumed to be perturbed by additive 
Gaussian noise, e, so that 





(Kendall and Stuart, 1967), 



Y = F(9) + e 



( 2 ) 



where e_ is characterized by an n x n moment matrix M; 



M = E { [ Y - F ] [Y - F^]' ) = E(ee' ) 
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The multidimensional Gaussian distribution of the error can 
be expressed as 

P (^) = C 2 exp ^ -1/2 J- 

or 

P(Y - F) = C 2 expj^-l/2[Y - “ f| j 

where C 2 is a positive constant (see, for example, p. 151 
Davenport and Root, 1958) . 

From equation (2) it can be seen that the probability 
density of Y given F is equal to the probability density of 
the error; that is 

P(Y(P) = P(e) = P(Y - F) 

and since F is a function of 9 this density is also the proba- 
bility density of Y given 0. The likelihood function, there- 
fore, is equal to the probability density function of the 
error, namely a multidimensional normal distribution, that 
is, 

L(Y;9) = P(Y|e) = P(Y - F (9) ) . 

Thus 

L(Y;0) = C2 exp^ -1/2 [x - fJ' [y “ Z. J } (3) 

The likelihood function is a maximum when the magnitude of 
the exponent, ||y - X J [X ~ Zj • ^ minimum. 

When the noise ^ has zero mean and the error in each 
measurement is statistically independent of the errors in the 
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other measurements the moment matrix M is diagonal 

.2 

M = E(ee ' ) = 



^ (j-2 0 



0 *CT2 

n 



and M 



-1 



VC5 . 



0 



•i/cr 

n 



where CX^ denotes the variance of the error in the measure- 

ment. The moment matrix, M, can be written in terms of weight- 
ing factors Wj^ and the standard error O' . The ma trix of 
weighting factors is denoted by V 



and 



where w. 

1 






0 





-1 

V 




0 




crVcr. 

1 



2 



(4) 



The maximum likelihood estimate in this case is the weighted 
least squares estimate since the magnitude of the exponent 

in equation (3) is proportional to the weighted sum of squared 

« 

errors : 
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n 

1 y 

1/2 [y - F(e)]* M-1[y - F(0)J = ^ -(Yi - fi(0))2 
When the variances of the errors are equal then the weighting 
factors, are equal to one, and the maximum likelihood 

estimate is the familiar least squares estimate which minimizes 
the sum of squared errors SS(^); 

n 2 

SS(0) = Z (Yi - fi(©)) . 
i=l 



Assume that a rather accurate estimate of ^ has been 
obtained and that jF(^) can be approximated by a first order 
Taylor series 



F(0) = F{0) 

dE. 

where X = 

a© 

element ^f • 



+ X A© , 



^ is an n X p matrix of derivatives with ij 
© 

and the errors A © = © - 6 are assumed small. 



When this expression for JF(^) is substituted in equation (3), 
and when the ^ priori probability density P (^) can be consider- 
ed to be a constant over the region of interest, the probability 
density of £ given Y can be shown to be a multidimensional 
normal distribution with moment matrix T : 



P(©|y) = C 3 exp[ J 



(5) 
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—1 J -2 

where J_ = X’M X = CT 

The moment matrix _Q is known as the variance-covariance 
matrix since it consists of the variances (diagonal terms) and 
covariances (off-diagonal terms) of the estimates. it is also 
known as Fisher's information matrix (Van Trees, 1968). Specif- 
ically, under the conditions of uncorrelated, zero mean Gaussian 

noise, the variance of the estimates can be expressed as 

1 -1 

11 >4 JL, A 

This estimate is approximately equal to the asymptotic 



n 


^ f 


L 


w. 

9 0./ 1 


. i=l 


1 / 



variance 1 /r (9) is close to 0. 



4) Application of Parameter Estimation Techniques to the Survey 
Problem 

The pulse-Doppler navigation system makes measurements 
which are characterized by a set of unknown parameters. The 
unknown parameters are the geometrical coordinates of the 
beacons and the ship positions at which the measurements are 
taken. The measurements consist of slant ranges and slant 
range differences. 

An orthogonal XYZ coordinate system will be used in which 
the beacon coordinates are (0 , 0, Z^^) , . (X 2 , 0, Z 2 ) , and {X^,Y 2 
as shown in Figure 2. (Only the 3 beacon survey problem will 
be considered here; the method can easily be extended to more 
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(XSi,YSi ,ZS) 
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Figure 2. Survey Coordinate 



2 







I 



i 

I 



I 






than 3 beacons) . The depth of the ship hydrophone is assumed 
to be a known constant, ZS. Simultaneous measurements of slant 
ranges and slant range differences are taken at m survey points 
(XSi,YS^,ZS) , 

The unknown parameters, that is, the beacon and ship posi- 
tion coordinates, can be expressed as a column vector £, with 
transpose 

» © 2 ' * * * ' ^ 

where p = 2m + 6. 

The "noiseless measurements" are the exact geometrical slant 
ranges and slant range differences. These can be expressed as 
= SR^j(9) = [(XS^-Xj)^ + (YSj_-Yj)2 + (ZS-Zj)^ ] 

where ^ ^2 ^ i — l,2,...,m 

j = 1,2,3 

Similarly, the slant range differences between two survey points 

i and ii and the beacon can be expressed as 

2 2 2 1/2 
DSRij(e) = [(XSj_i-Xj) + (YS^j^-Yj) + (ZS-Zj) ] 

2 2 21/2 
- [ (XSj_-Xj) + (YS--Yj)^ + (ZS-Zj) ] 

These functions can be expressed as a column vector ^(^) 

with transpose, 

F’ (£) = [ ^^3' ^^21' • • • '^^^ll'^^^12 '^^^3' ' * * ^ * 

I 

For m survey points F(^) has dimensions 3m + 3 (m-1) = 3(2m-l). 

If F(^) could be measured without errors then some minimum 
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number of measurements would be sufficient' to determine ^ 
exactly. In this case the slant range differences are linear 
combinations of the slant ranges and the minimum data case, as 
shown by Vanderkulk (1961), requires 18 slant range measurements: 

SRj^j(e), i = 1,2,. ..,6, j = 1,2,3 
or equivalently, 3 slant range measurements and 15 slant range 
differences 

S%j (9) , j = 1,2,3 

DSR^j (9) , i = 1,2,. ..,5, j = 1,2,3 . 

If only slant range differences can be measured then 24 measure- 
ments, 

DSR^j(9), i = 1,2, ...,8, j = 1,2,3 
are sufficient to determine ^ exactly provided one has reasonable 
initial estimates of the parameters. This condition is a result 
of the linear approximation used to solve the non-linear equations 
The observations Y differ from due to the presence 

of noise. Slant range observations from the transponder system 
can be expressed as 

Sij = SRij(£) + €ij 

where errors, ^ij' assumed to be independent identically 

distributed normal random variables with variance • Obser- 

vations of slant range differences from the Doppler system can 
be expressed as 

DS. . = DSR. . (0) + S. . 

ID ID - ID 



2 5- 




% 






where the errors 



, ^ i j ' assumed to be independent identi- 

2 

cally distributed normal random variables with variance Cf ^ . 

The errors for the complete system can be expressed as 

f 



e = 



L(^J 



= [X - F(t)] * 



If the errors ^ and 5 are independent of each other^ than ^ 
is characterized by the moment matrix 



M = VcX = 



cTiV^^ 0 



0 






^i" = 

i = 1,2, 
^m+i = 



, 3m 



i = 1, 2 , . . . , 3 (m-1) 

5) Gauss-Newton Method of Non-linear Least Squares Estimation 
We have seen that the maximum likelihood estimator under 
conditions of additive, uncorrelated, Gaussian noise with zero 
mean is equivalent to the weighted least squares estimator and 
that this estimator is consistent and' efficient. The weighted 



least squares estimate is the ^ which minimizes the residual 
sum of squares 

SS(0) = e'V"^e = [x - Z(9) ] [X “ Z(Q) . (7) 

In general the estimate ^ is the solution of the p normal 



equations 



9 SS(9) 

= o i = 1,2, ...,p 

9 ©i 

« 

assuming that the solution is interior to the parameter space. 
When these equations are linear they can be solved directly. 
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when the normal equations are non-linear as they are in this case 
an iterative technique must be used. A number of techniques 
exist; for example, quasi-Newton, conjugate gradient, etc. A 
brief review of numerical techniques for fitting non-linear 
models is given by Chambers (1973) . 



The method used to find the estimate £ is based on a modi- 
fied Gauss-Newton procedure. The Gauss-Newton procedure consists 
of linearly approximating the function F(_0) and applying standard 
linear least squares techniques. 

One linear approximation of F(^) about an initial estimate 
^ IS a simple Taylor series 

F(9) = F(^) + X Ae 

9 F 

where X = a is an n x p matrix of derivatives with 

- 90 ^ 



elements x^j 








and 




/N 

9-9 is a column vector 
- -0 



of length p. Choosing the new variable Z = Y - F(^) we have 
the linear model 

Z = X A9 + e. 

The linear approximation to the residual sum of squares is 
SS(9) = e'V“^e = j^Z-X A 9 J'V“^ J" Z-X A9 j= SS(^) , 

The Gauss-Newton step A 9 is obtained by linear regress- 
ion techniques on this linear residual sum of* squares 
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I 



(Draper and Smith, 1966): 

A9 = (x'm“^x)"^ 

A /A. 

A new estimate ^ +A9 can be formed and the process 

repeated, until a given convergence criteria is reached. 

6) Modifications of the Basic Gauss-Newton Procedure 

Modifications of the basic Gauss-Newton procedure can be 
made to insure that the step A 9 , will reduce the residual sum 
of squares and avoid singularity problems in the 
matrix. A complete discussion of the modifications used/ 
based on stepwise regression techniques, is contained in 
Jennrich and Sampson (1967). A summary of that discussion is 
presented here. 

One important modification of the basic Gauss-Newton 
procedure is the use of partial steps, A 9 , in place of the 
full step A 9 when the full step results in an increase in the 

A 

residual sum of squares SS (£) . This frequently occurs when 
the linear approximations upon- which the method is based fail. 
The proportion to be used is obtained by trying the arbitrary 
sequence of values 1,1/2, 1/4, 1/8,... until the residual sum 
of squares is reduced. It can be shown that a sufficiently 
small step in the direction of A 9 will always reduce the 

A 

residual sum of squares unless ^ is already a minimizing value. 
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Another modification of the basic procedure is dictated 
by the possibility that the matrix is, or is nearly, 

singular. Step-wise regression techniques are used where the 
parameter selected at a given step is the one which makes the 
greatest reduction in the residual sum of squares. If there 
is a singularity problem then only a subset of the parameters 
is used. This also provides a convenient way to handle bound- 
ary restrictions on the parameter values. A parameter is 
modified only if the new parameter value is within or on the 
parameter boundary. 

7) Iterative Technique 

In summary the procedure used consists of the following 
steps : 

1) The value of the function F(^) and its derivatives X 

/V 

are calculated using an initial estimate . 

2) The X'X~^2i i^atrix is formed and the residual sum of 
squares SS (0Q ) is computed. 

3) The X'Y."'^2i matrix is inverted in a step-wise manner 

and the step size A 9 = (X'Y ^X) ^ formed. 

A 

4) A new estimate ^ + A 9 is formed and the new 

residual sum of squares SS(^j^) computed. If SS(^j^) is larger 
than SS(0 q) then the step size is halved until SS(9j^) - SS(0q) 
or until 7^ = (1/2)^*^. 



- 29 - 



5) This completes an iteration. The process is repeated 

until the solution converges. Convergence is determined by 

the relative change in the residual sum of squares, i.e. the 

solution has converged when the relative change is less than a 

specified criterion, C, 

|SS(9n+l) - SS(9j^)| 

■ = C'< C . 

SS(9„+l) 

6) In practice we require that 4 successive iterations 
result in C'< C in order to accept the hypothesis of conver- 
gence. In some cases the algorithm may fail to converge to a 
minimizing solution (e.g. when it converges to a local, rather 
than a global, minimum) . A plot of the parameter estimates 
and examination of the residuals will usually provide some 
indication of this. 

8) Error Estimates 

We have seen that maximum likelihood estimates are 
asymptotically normally distributed with known variances and 
that these variances are the diagonal , elements of the variance- 
covariance matrix _T. For a process perturbed by uncorrelated 
Gaussian noise with zero mean this matrix is given by the 
expression 

r = (X' V"^X)“^C7'2 where = 



wi 0 

Wo 

0 ' * 



n 
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The (X'Y“^X) niatrix is called the correlation matrix and its 
inverse will be denoted by A. 

til. 

The standard error in estimating the i parameter, O 

0 . 

1 

is the square root of the i^^ diagonal element of the variance- 



covariance matrix 



T: 



(X. = o ' VaT 



( 8 ) 



th 

where a^ is the i^ diagonal element of the A matrix 



ai = 



n f 






2 


t — ■■ ' 


3e . 
1 


A 

9 / 





-1 



(9) 



The terms on the right hand side of equation (8) are indepen- 
dent. The standard error, O' , depends upon the accuracy with 
which the slant ranges and slant range differences can be 
measured. An estimate of the standard error can be obtained 
from the residual sum of squares 

Cf = SS(0) 

n-p 

A 

where SS (9) is defined in equation (7). 

The other term, V a . , is sometimes called the error 

1 

magnification factor and is a function of the survey geometry. 
The magnitude of the error magnification can be obtained by 
evaluating equation (9) . 
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The approximate errors in estimating the beacon coordin- 
ates 

( 7 ^ — aj ' ^ ~ 1 , 2,. ..,6 

^i 

can be reduced in two ways: 1) improvements in the measurement 
accuracy and 2) a suitable selection of the survey coordinates. 
In the next section, we will consider in detail how this may 
be accomplished. 
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Ill MINIMIZATION OF ERRORS 



1) Measurement errors (pulse system) 

Slant ranges and slant range differences are measured by 
two different systems for which the magnitude and sources of 
error will in general be different. A discussion of these 
errors requires a complete description of the system operations. 

The pulse system is designed to measure round trip travel 
time of an acoustic pulse between a ship-mounted or ship- 
suspended transducer and near-bottom transponders. When an 
acoustic pulse is transmitted from the transducer, digital 
counters in a shipboard receiver begin measuring elapsed time. 

The transponders receive and recognize the transmitted pulse 
and generate a "reply" pulse at a specific frequency. When 
a return pulse is detected by the shipboard receiver at a 
given frequency the corresponding counter is stopped and the 
elapsed time is displayed and transferred to a digital com- 
puter for processing. 

There are fixed time delays associated with this process; 
e.g. signal recognition time, transponder turn around time 
and signal processing time. Accurate calibration of the system 
can reduce the uncertainties in each of these delay times to 
approximately ±.5 milliseconds (msec). For most applications 
these timing errors can collectively be considered to be normally 
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distributed with zero mean and standard error of 1 msec. Two 



i of these timing errors, signal recognition time at the trans- 
ponder and at the ship are functions of the signal-to-noise 
ratio and the pulse bandwidth. For example, the errors in 
recognition time increase to approximately .8 to 1 msec for a 
signal-to-noise ratio at the pulse receiver of 27 dB. This 
signal-to-noise is typical of a range of 10 km in sea state 3 
for the system presently in use (source level; 189 dB, re: 

I yu Pa pulse length: 10 msec). For ranges greater than about 

II km (for the present system configuration) signal recognition 
is erratic and the errors can be assumed infinite. For fixed 
bandwidth, timing errors can be reduced by increasing signal 
power or they may be acounted for by somd appropriate weighting 
of the data. 

Another source of error in the measurement of travel times 
is ship motion. The error in round trip travel time due to 
horizontal ship motion is given approximately by the expression 

€ = TT^V^COS e/c 

where TT = measured round trip travel time 

0 = angle from the vertical between the transponder and 
ship transducer 

Vjj = horizontal component of ship speed away from the 
transponder 
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c = average sound speed 

For ship speeds of 1 knot and measured travel times of 6 to 
12 seconds (typical in deep ocean applications) the round trip 
travel times will be in error by approximately 2 to 4 msec. 

This error can be reduced by taking measurements while stopped 
or by reprocessing corrected data once' the approximate ship- 
positions and speeds are known. Other ship motion such as 
heave, pitch and roll cannot easily be accounted for but when 
a large number of measurements are used their effect can be 
considered as an uncertainty in transducer depth equivalent to 
a timing error of approximately 1 msec. 

A summary of these errors is present in Table 1. To achieve 
these errors in measured travel time of a few milliseconds 
it is essential to use an accurately calibrated system and 
make the measurements while stopped. These errors can be 
assumed to have zero mean for a large number of measurements. 

To obtain geometric slant range from corrected travel 
time the velocity of sound must be known. Although the sound 
velocity has spatial and temporal variations, horizontal and 
temporal variations can be assumed negligible. Sound velocity 
variation over depth, however, is commonly on the order of 
50 meters/sec for 5 km deep water. The structure of this 
variation, the sound velocity profile, can be determined from 
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Table 1. Sources of errors in measurement of Travel Time 



Source of Error 


(f (msec) 


O’ (meters) 


Recognition Time 


1.0 


1.5 


Miscellaneous Timing Errors 


1.0 


1.5 


Horizontal Motion < 1/2 knot) 


1.0 


1.5 


Other Motion 


1.0 


1.5 


Total = 


2.0 


3.0 
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historical data or direct measurements. Historical data can 
be expected to be in error by 2.5 m/sec or more, whereas 
present measurement techniques have an accuracy of about 
.25 m/sec. For a precision navigation system the use of direct 
measurements is strongly preferred. When the sound velocity 
profile is known the travel time of an acoustic signal between 
two points can be determined and an effective sound speed, V^, 
can be calculated. 

Vg = SR/TT 

When the travel time is measured using a pulse navigation 
system the effective sound speed is not known since it depends 
upon the positions of the transponders and ship transducer. 

The slant range, however, can be determined approximately by 
using a calculated range R, equal to the corrected one-way 
travel time TT^ multiplied by the arithmetic mean sound 
velocity c’: 



SR R = CTT 

c 



1 

where "c = 

Z . - ZS 
1 




c (z ) dz 



The depth of the transducer ZS is usually known and errors of 
less than 50 meters in the estimate of transponder depth, Z^, 
have a negligible effect on c’. The arithmetic mean sound 
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velocity can therefore be considered a known constant. 

The ratio of slant range, SR, to calculated range R = "cTT^, 
can be expressed as a function K(Ve) : 

K(ve) = Sr/r = ve/c 

For a given sound velocity profile and depths ZS and Z^ the 

effective sound speed, Ve, is a monotonically increasing 

function of the geometric angle between the transponder and 

the survey point measured from the vertical (Vass, 1966) (see 

Fig. 3). Since the calculated range, R, is also an increasing 

function of the geometric angle, K(Ve) can be expressed as a 

function of R and approximated by a second order polynomial 

2 

K(Ve) = K(R) ~ aQ + a^^R + a 2 R . 

The coefficients, a^, are found by calculating the range R 
and the ratio SR/r for various angles and applying a second 
order curve fitting subroutine. Travel times are computed 
using standard ray acoustic techniques (Officer, 1958) . A 
description of program ARCOR which performs these computations 
is contained in Appendix A. The error in using the relation- 
ship SR = R is usually less than 1.5 meters but can be as 

great as 3 or 4 meters for angles greater than 80°. When the 

2 

relationship SR = ^ used the errors in 

slant range estimates are less than .1 meters. When the sound 
velocity profile is accurately known, errors caused by acoustic 
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refraction are negligible compared to timing and transducer 
motion errors. 

2) Measurement Errors (Doppler system) 

The Doppler system in its present configuration (Fig. 4) 
consists of three bottom-moored CW beacons each separated by 
50 Hz, a phase detector for each beacon and a real time com- 
puter processor. Each beacon contains a crystal oscillator 
with a frequency stability of 10"^ Hz and has a nominal source 
level of 166 dB (re ±1 Pascal (|) 1 m) . The beacon signal after 
being processed by its phase detector is recorded on analog mag- 
netic tape and is fed to a digital computer for real time pro- 
cessing. The phase detector output consists of two base-band 
outputs one proportional to the sine of the input phase angle, 
the other proportional to the cosine of the phase angle. Phase 
variations due to a finite receiver velocity v(t), can be written 

as 0 {t) = 2 TTf tv(t)/c. If the input signal is represented by 
B 

s (t) = A (t) exp [ -j [ 2 TTfgt - 0(t) J j 
the phase detector outputs are given by 

Sj^(t) = A(t) cos 0’(t) 
and S 2 (t) = A(t) sin 0 (t) . 

The instantaneous phase of the beacon signal can be determined 
from the ratio 

s (t)/s (t) = tan 0(t) ■ 

2 1 . 
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Figure 4. W.H.O.I. Doppler Receiver 



In practice the quadrant into which the phase angle falls is 
determined simply by examining the signs of Sj^(t) and S 2 (t) . 
Phase differences between phase detector samples are computed 
and accumulated in quanta of A /4 3 cm at 13 kHz. The 

accumulated phase difference after T seconds is of the form 



where f is the sampling frequency of the phase detector out- 
s 

puts (10 Hz) (Porter ^ » 1973). 

The distance, DSR, travelled toward or away from a beacon 
in time T is the accumulated phase 0 ^ multiplied by the wave- 
length Xof the signal: 

DSR = A 

The error in measurement of accumulated phase depends upon the 
number of samples N = Tf_ and the signal detection probabili- 

o 

ties of the receiver. Since the samples can be considered to 
be independent measurements the variance of the accumulated 
phase is the sum of the variance of the individual samples 




N 



z 

n=l 



2 




2 

N(X 

0n 



ana f(p^) 

* tin 

where Pe is the probability of error in the n sample, 
probability of error is a function of the signal-to-noise 



The 
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ratio at the phase detector (see Appendix B) . Table 2 shows 
some representative values of Pe in a . 1 sec interval for 
various signa 1-to-noise ratios encountered in practice. The 



both number of quadrants and meters. A quarter wavelength of 
3 cm is assumed. Nominal horizontal range is calculated for 
a source level of 166 dB (re: lyc^Pa) and a noise level of 
44 dB/Hz (re; 1 j^Pa) assuming transmission losses are due only 
to spherical spreading (20 log r) and attenuation (1 dB/km) . 
Sample calculations are given in Appendix C. 

Random phase fluctuations can be caused by multipath 
interference and forward scattering. The root mean square 
phase fluctuations from these sources has been estimated to 
be about .13 quadrants (Porter £t , 1973). When compared 
with the errors due to ambient sea noise shown in Table 2 



Slant range difference is obtained from accumulated phase 
change by multiplication by a scaling factor: the wavelength 



primarily by errors in sound velocity (the instability of the 
beacon frequency is negligible) . When the arithmetic mean 



associated errors and (T = 1 hr) are also given in 



they have little effect (about 




and 



can be considered negligible. 



of the signal. Errors in wavelength X = c/fg are caused 
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TABLE 2. Doppler Errors vs Signal-to-Noise Ratio 



10 log S/N 


P 

e 




0^ (1 hr) 


X 


(dB) 




(quad) 


(quad) 


(meters) 


(km) 


33 


.023 

1 


.217 

1 


41 


• 

to 


1 

0 


30 


1 

.031 1 

i 


1 .253 




48 


1 


3.7 


27 

1 


.042 ! 


.293 


56 


1.7 ^ 

i 1 


5.5 


I 

24 


.055 i 

i 


.341 

1 


65 


1 i 

1.9 


7.5 


i 21 


.072 


.393 


75 


2.2 


9.2 


18 


.095 


.456 


87 


2.6 

1 


11.1 
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sound velocity of an accurately measured sound velocity pro- 
file is used, the standard error in wavelength will usually 
be less than 1 meter. This can be reduced further by using 
the average sound, 

S = length of travel path/travel time, 
once the approximate survey positions are known. 

In summary, when operating above the pulse system thresh- 
old/ slant ranges can be measured by the pulse system with an 
accuracy of 3 to 5 meters. Under the same conditions slant 
range differences can be measured with an accuracy of 1.2 to 
2.5 meters for measuring intervals of 1 hour. 

3) Selection of Survey Coordinates 

The error in estimating beacon coordinates also depends 
on the error magnification terms of the beacon coordinates, 

= 4 ^, i = 1,2, 6 

where the a^ are the i^ diagonal elements of the A = (X V X) 
matrix. These terms are dependent upon the number and relative 
locations of the survey points. In general the error magnifi- 
cations will vary inversely with the degrees of freedom which 
increase with the number of survey points above a certain 
minimum. For example the slant range difference equations 
require at least 9 survey points for a solution to exist, and 
the degrees of freedom for this set of equations is equal to 
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the number of additional survey points. Comparison of survey 
techniques should be based on the same number of degrees of 
freedom. For example a pulse survey of 7 points (1 degree of 
freedom) should be compared to a Doppler survey of 10 points, 
and a pulse survey of 10 points should be compared to a pulse- 
Doppler survey of 10 points. 

A measure of "good" survey locations is the trace of the 
beacon coordinate covariance sub-matrix, that is, the sura of 
the squares of the error magnification factors of the beacon 
coordinates : 



TRACE = 

i=l 



Mi 



6 

= E 

i=l 



The trace is a function of the beacon and survey location 
coordinates, i.e. the parameter vector and will be denoted 
by TR (^) . For a given number of survey locations there will be 
a set or sets of survey locations for which the trace is a mini- 
mum. A typical deep ocean deployment of the pulse-Doppler 
navigation system will be discussed for illustration. 

The most commonly used 3-beacon acoustic navigation net 
is ideally an equilateral triangle. The baseline length, i.e. 
the length of the sides of the triangle, are typically on the 
order of the beacon depths: about 5 km. For this configura- 

4 

tion of the net the covariance matrix was evaluated for several 
sets of survey locations, each consisting of 10 survey points. 
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A few examples are given in Figures 5A through 5E. The value 
of the trace of the beacon covariance matrix for each arrange- 
ment is listed in Table 3. The trace for the corresponding 
survey with only pulse data is given for comparison. A more 
rigorous technique for determining "optimum" survey locations 
was also used. Powell's algorithm to find a (local) minimum 
of a function of several variables was applied to a function 
G(9) = TR{9) + C (^) , where C{9_) is an arbitrary function added 
to the trace to constrain the survey locations within the 
maximum range of the system (Powell, 1965) . No significant 
improvement to the starting value of G(9) 3.6 was realized 

although several survey geometries were tried. 

4) Computer Simulation 

A computer program called PDSRV was developed to solve 
the survey problem using the parameter estimation techniques 
discussed. It is a modification of a Biomedical Computer Pro- 
gram (BMD 07R) (see Dixon, 1973) . The program requires pulse 
data from at least six survey points and/or Doppler data from 
at least ten survey points. In addition to the data, the 
user must input the average sound velocity and the acoustic 
refraction coefficients found by program ARCOR. The beacon 
frequencies, the depth (s) of the transducer and initial 
estimates of the beacon positions must also be provided. 
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TABLE 3. 



Comparison of the Trace of the Beacon Covariance Matrix 



for Pulse/Doppler and Pulse Surveys 



Figure # 


Trace 

Pulse & Doppler 


Trace 

Pulse Only 


5A 


3.6 


6.8 


5B 


4.0 


8.3 


5C 


4.3 


10.4 

f 


5D 


17.0 . 


58.2 


5E 


138 


325 
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Maximum and minimum values for the beacon coordinates are 
constrained to provide more rapid convergence. The constraints 
are set as follows: 



X . 


= X. 


±2000 


m 


3 


3 






Y. 


= Y. 


±2000 


m 


3 


3 






Zn 


= Z. 


+ 1000 


m 


3 


3 







where (X j , Y ^ , Zj) are the initial estimates of the beacon 

' coordinates. Other parameter bounds may be specified if desired. 

Although it is unlikely, even with the most rudimentary primary 
I navigation system, that initial estimate errors will exceed 
I those specified above. 

The program consists of three main sections: (1) initiali- 

zation and data input; (2) minimization of weighted squared 
, error; (3) determination of the • approximate error covariance 

I 

I 

I matrix of the beacon positions. For convenience the program 
I is divided into 7 subroutines as illustrated in Figure 6. 

1 A description of each subroutine is given in Table 4. 

Two additional subroutines were used to simulate survey pro- 
blems. Subroutine READ calculates data for a given set of survey 
I locations and beacon coordinates. Subroutine NORMA generates 

I 

Gaussian distributed numbers with a specified mean and standard 
^ deviations. These are used to simulate data with errors. 
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TABLE 4. Description of Program PDSRV 



Name of 
Subroutine 



MAIN 



DINPT 



MINIZ 



XPRMX 



Description 



Reads the following user inputs: 

1. Transducer depths 

2 . Convergence criteria 

3. Bounds for beacon parameters 

4. Beacon coordinate estimates 

5. Survey position estimates (optional) 

6. Bounds of Survey position (optional) 

1. Reads input data from specified device (s) 

2. Calculates slant ranges (SR) from travel 
times 

3. Calculates changes in slant range (DSR) 
from accumulated phase dhange 

1. Computes the mean squared error (EE) 
from the weighted sum of squared errors: 
EE = 1 /DOf) SS(^) 

2. Checks that EE is less than previous 
value; if not the step size is halved 

3 . Checks convergence criteria 

1. Computes the elements of the X V ^ X 
matrix 

2. Calculates the weighted sum of square 



errors : 




I 




I 



1 



TABLE 4. 



Continued 



Name of j 

Subroutine 



i 



I 



FUN 



STEP , 

I 

I 



1 OUTPT 




Description 



1. Calculates values of the slant range 
and slant range differences (fj^^) using 
latest estimates of beacon and ship 
positions 

2. Calculates the elements of the X matrix, 
i.e. the derivatives of the slant ranges 
and slant range differences: 

1 

X . . = — ^ 

d 0j 1 0 



1. Inverts the X V^~X matrix in a step-wise 
fashion 

2. Computes the step sizes A 6 

Computes and outputs the following 

1. Final. estimates 

2. Correlation matrix for Beacon Positions 

3. Root mean square error (standard error, cr' ) 

4. Error magnification terms of the 
covariance matrix 

5. Residuals and standard errors of each 
measurement (optional) 
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Typical output of a simulated survey problem is shown in 
Tables 5a through 5d. Table 5a lists the beacon and survey 
coordinates used by subroutine READ and the errors in data 
generated by subroutine NORMA. Table 5b lists the simulated 
data which is input to program PDSRV. The final estimates, 
correlation matrix, error magnification factors and standard 
errors of the beacon coordinates are also shown in Table 5b. 
Similar information for the survey locations is listed as 
illustrated in Table 5c. The residuals, 

SR(observed) — SR (estimated) and 
DSR(observed) — DSR (estimated) 
can be compared with the estimates of standard error of each 
measurement. This output, labeled SURVEY' ERRORS, is optional 
and is shown in Table 5c. 

After each iteration of the minimization process the mean 
squared error and latest estimates of the beacon positions 
can be output as shown in Table 5d. The number of the iteration 
and the number of times the step size was halved is also 
included. 

Survey problems with measurement errors of more than 750 
meters have been simulated. The minimization process converged 
in every case and the beacon coordinate estimates were nearly 
normally distributed about the true values with the standard 
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U1 1 . 800 



Errors in the slant range 
neasurements to Deacons 1,2, 
and 3 ) Tiaae at each survey point. 



ERROR DOPPLER DOTH 



16^027 

2. 822 
5. 345 
-7, 257 
13.01? 
. 165 

C*. ] J 


- IK 181 - 

-3. 546 
-IK 357 
-12. 374 
3.525 
- 3 . 702 
S ^^47 


K 064 
1 1 . 563 
- , 264 
- 1 5 . G. 1 6 
-14, fi93 
4 . 3G".6 


Errors in the slant range 
difference measurements 


-4.913 


- 7 . 530 


-8. 272 


tjetween survey points 
142, 2 & 3 , 3&A, etc. 


- -7.140 


-14. 132 


,751 




MERH= . 600 


STD. DEU.= 10, 


00 UT= 1 . 000 





Table 5a. List of beacon and survey coordinates used for survey 
simulation. Random normal errors generated by Subroutine NORl-.'A. 
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— PULSE -SV3TEM- Dm R Slant Ranges to Beacons 

HR TIHECHcMiS) 1 2 

— SLR 1 F.?85t-1-} 

SLR 2 142*«:M2-*:142 7695.95 7699. S3 

— SLR --3 213*213+213 7674. 9S - 10675.67 

SLR 4 2S4+264+2S4 5687.65 7601.89 

—SLR — 5 —--355+355*355 7655.63 9978.68 

SLR 6 426+426+426 7500. 98 7491.57 

— SLR 7 497+4 97*497 998 1-. 6 1 7059.-73 

SLR 0 568+56S+56S 7608.33 ' 5105.27 

— -SLR 9- — ■ 639*C.39-*639 - 16665.83 7666,11 

SLR 10 710+710+710 11465.76 11453.79 



5112. 7 
-7457.5 
7765.6 
1 I 167. ^4 
165 ^ 8 . 69 
4 1 1-4 7-. 68 — 
7763. 97 
.7444. 97 



7565. 75 



— 




ii n n T n 


Slant *-Range 


Differences" vdth respect to Beacons 




HR 


TIHE<H+M+S> 


1 


2 


3 


- DSLF‘ 


1 


- 71+ 71 * 71 





.00 


• 06 


DSLR 


2 


1 42+142+142 


I8‘98.50 


1926. 51 


“651 . 05 


C'SLR 


3 - - 


-213+213+213 


--19. 28 


2970. 29 


2381 . 9J 


DSLR 


4 


284+284*284 


“2574. 74 


“3074, 02 


317, 3? 




r.“ 


'“i S *^ *^. *^ * *^ S - - - 






*T*?r ft r , r , 




V' 




J ^ H Cl 


— 


-o >> t • . O 


PSLR 


6 


426i=4 26'i426 


460. 36 


“2433.47 


“585. 70 


68LR 


7 


4 9 7 4 9 7 + 4 9 7 


248?. 16 


“451 . 04 


575.91 


OSLR 


8 


568-»^5e.8i'56S 


“23?8. 34 


“1949. 65 


“3376.54 


- f'SLR 


Q 


639»t«639^-639 


3657. 71 


2576. 56 


“325.91 


DSLR: 


1 6 


7 IfH:? 16=4^7 16 


772. 22 


3761 . 1 6 


97. 27 










.■M • - 







r IHRL 



tST IHRTCS 




Y 

. 66 



5668.69 .60 

25 i 4325^-9^ 



Z 

4995,86 
5663. 22 
16.- 



( 2 ) 



cor<Ft:LnT JOH 


nniRiK For< 
X2 

1 . 0000 
, 4 Q H 6 


BETnenw 

.X3 

. 4 986 
1 , Anmi 


POSIT I OMS 
V3 
1951 

*^04 4 2 


21 

“. 1030 
0S04 


22 

“. 1242 
“ . 6872 


23 

. 1 263 
.6970 




1951 


“.0442 


1 . 0000 


. 1299 


. 120? 


. ti339 


- ... - 


1830 


“, 0SG4 


1299 - 


1.U000 -- 


. 8109 


.1588 




“ . 124 2 


0072 


. 1 207 


.0183 


1 . 0000 


.1313 





- . 1263 - 


6970- 


. 0339 


-.1590 


.1313 


1 . 0066 








(3) 








LKF,- HRCi, 


1^0947 


1 . 1 


, 7 1 4 ' 




. r.3i2 


5850 


STD, LkR. 


10.9884 


1 1 . 3696 


6. 7401 


6. 0128 


5.934 8. 


5. 8725 


_ TKOCE 


_ 


3. yp3 










STRHOnfv'D ER. 


ROR IN DOPPLER DmO 


- 10.038 











. ... 





(A) 


- - 


— 





Table 5^>. (1) Simulated Pulse f^Doppler measurements. (2) Final estimates of 
beacon positions, (3) Correlation matrix for beacon positions. 
(4) Error magnification factors and estimated standard errors for 
each beacon position. 
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SURVEY. POSITIOHS _ 



_HR X ERR.HftG STD. ERR V ERR.MRG STD. ERR Z 

1 1 .-53 14 92.-20 1 .-37 J3. 7lC«7C . 60- 

2 2495.33 1.57 15.74754 5292.56 .66 6.66614 .66 

___5 *^2939.46 - 1.45-14.55447 4997.56 1.36 13.66243 .60 

4 -676.73 1.37 13.72696 -711.37 1.47 14.72249 . 6Gi 

_5 -2463.60- 1.96 19.7Gf095 4317.24 1.56 15.65258 . 6C; 

6 2522.26 1.55 15.533.95 -4997.62 .68 6.30346 .60 

7 7502V79 1-.5.(3 1 5. 8 14T-5 1 =4299^39 l-,,-5f 1 5 , 56735— r OG'— 

8 5697.97 1.05 10.55979 -698.21 1.46 14.67303 .GO 

— 9 7995.33 1.25 12.57312 - 4996.21 - 1.36 13.62897 .00 

16 2531.10 2.68 26.86143 9991.16 .62 6.19808 .60 



( 1 ) 

— SUWEY-F.KFGK5 



— 


residual: 


5 <0BS ~ 


EST> 


- STD. ERROR 


OF ESTIMATES 1 


NORM. 5E. 


PT 




BEfiCOH 






BE room 






HR 


- \ 


^ _ O 


3- 


“* 


2 


3 




1 


C .• £ Q O 




« 4 , 3 


7 -282 


*7 — » 


7_, 54^; 


\ , 0643 




- — - ” v> * O I? ■ ' 














2 


2. DS7 


-5.346 


16. 997 


6.713 


6.719 


5. 768 


1 . 7868 


.. 5 


2.254 


. 246 


-18. 566 


6 . 63? 


6. 273 


6.467 


1 . 0766 


4 


-3. 725 


-4 . 394 


-18.124 


5. 924 


6 . 680 


6. 672 


1 . 8945 


-- 5 


- 8. 223 


-3. 855 


9.355 


6.170 


6.676 - 


6. 185 


1 . 2989 


6 


-2. 887 


-4 . 2S3 


16. 977 


6. 394 


6. 356 


5. 984 


1 . 7565 


— 


“‘Sr 121 


— V .~6 3 o 


2, 545 


6-. 67 6 


- 6 V 1 1 3 


6. 6 72 


, 9034 


8 


-1 . t.9S 


6 , 673 


-8. 429 


6.691 


5, 895 


6. ( 25 


1 . 084 1 


9 


- -4. 3S1 


-9 . 884 


-6.652 


6 . 368 


6, 752 


6. 4 15 


1 . 2751 


J6 


1 1 . 98U 


8. 668 


7. 521 


7. 54 5 


7. 570 


6. 1 38 


1 . 6523 




-1 . 579 


2.861 


3. 5 65 


7. 707 


7 , 675 


6, 94 7 


. 4532 


3 


1 . 83c. 


- . 961 


9.5)8 


7.418 


7. 109 


6. 792 


. 9? 1 4 


jt 


iT tl Ct O 


ri- gi 1 V. 7 l«l 




,1 i O 


7 ^2^7- 


7 ,365 


, .023 i 





— Of C'C* ->■ 















5 


-7. 892 


-8. 383 


-14.278 


5.169 


7.409 " 


6. 892 


1 . 8180 


6 


4 . 781 


-4.496 


-8. 238 


6.641 


7.310 


6. 327 


1 . 0450 


? 


3. 4 93 


-7. 898 


7.4 82 


7. 330 


6. 644 


6. 365 


1 . 132? 


8 


- 1 . 839 


4 . 449 


. 687 


7.419 


5.831 


6.951 


, 4772 


9 


“2. 488 


-.84 5 


-5. 327 


7. 274 


6. 710 


7. 4 06 


. 591c> 


— 18- 





1 

V 

1 


-9t:-4-34 


7-.-6LV9 


— 7-t4^4 


—7^ 093 


\ e650O 



( 2 ) 



Table-5c, — ( 1 )- Estimates or survey positions- (XjYj vd.th their associated 

error magnification factors and estimated stamdard errors. 

(2) Residuals of slant range measurements (1-10) and slant 

range difference measurements (2-10). 
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2 £» 7 . 2b" 



5 0 ' 1.36017E+06 

'3 H b 1“ 6 5 ^ 66 4 o 3 P's* 6 2" 



1 . 36778E+06 



8 0 1 . Bb873E-< GC 

— - •* ~ 

^ :» «:* . c» d 



f A c» • oo 



» c-> ^ H . d- 
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-trl4*4v5t? 



M 0 - - -S.03363E + G3 



12 6 6.21777E+G3 

4’b7v»-r-7 6 4-bGb^ bd 3V->9 r-l-c' 5 f^‘-4ri'r-4-* 



13 O 5.1G3mE + G3 - 

474J.96 35b5.G9 3625.33 5554.66 
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2624.^6 3664^?^ 
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-13 -^ -1 i . 156‘IC.EH 0 3 

4626.69 2337.39 



16 0 



A . 0671 5E+03 



•y 

oo 
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“» O C' V* • c» c» - 

46S6. 23 

A 'jtT tr '-i 


VI vr 


^VT JL. .“OZ 




try oo o.-. 


Cfc# 


oo OU . flCJ 


G<w - - 


—4686. ?9 “ 




* ~i o *? 'r 


2"ti 


H 5^ C* O . .*> O 


46 


— 5b bo. 5 9 — 


< 1 


ir .-> ,> r.' r •“. 


1 jt 


dO tl*. C-.i 


b ’ “ 


4c>3 1 ^ ^b' 


••: .< 
‘t n 


494J . 72 


n — 


-5201 .27- 


2i 


5679. 1 2 




c ~? ''fC ^ 




-- OOtflO. oo- - 


16 


5b BG. 66 


67 


— 5666 . 66 


4? 


5357, 6i 


1 Z 


5623 . i 9 



17 O 



1 . 01 663E -^0Z 



5606.91 


2512. 76. 


'1326. 17 


4S/SJ3. 6S/ 


5003. 6.5 


S.C1G3, 33 


J 0 0 

5603«69. 


1 . O07C.7E + O2 
2514 • iL. 


--•132S. S/2 


4SfSJ5.aa . . 


Su03. 22 


S603. 16 



Table Jd. Beacon coordinate ostljiiatoo (X2,X3>M3>ii1>Z2,Z3) at each iteration. 
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deviation approximately the same as that predicted. 

Single simulations with measurement errors of 5 to 800 

meters were made to determine the stability of the estimates 

A 

of standard error. O' , and the approximate error covariance 
matrix. Table 6 lists some examples showing the simulated 
measurement error, O ■ the estimated standard error, U , and the 
estimate of the TRACE of the beacon error covariance matrix. 

The stability of the trace indicates the linearity of the model 
in the region around 

The stability of the trace has another important implica- 
tion: suppose some "optimum" set of survey points have been 
determined for a particular navigation net; the actual survey 
points used might differ from the optimum by hundreds of meters 
due to navigation errors while conducting the survey; the 
resulting errors in the estimates of beacon coordinates, how- 
ever, would not be very different. The primary geometric factor 
which affects the errors of the estimates is the distribution 
of the survey points not their exact location. 
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TABLE 6 



cr 


A 

cr 


Trace 


5 


4.8 


3.96 


50 


47.8 


3.95 


500 


513 


4.76 


600 


640 


3.92 


700 


751 


3.98 


800 


820 


4.35 
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IV DISCUSSION 

Using the parameter estimation technique and the procedures 
described, the relative positions of acoustic navigation beacons 
can be determined with an accuracy of ±1 meter. 

A non-linear model is developed for a combined pulse-Doppler 
navigation system whose measurements of travel time and accumula- 
ted phase are perturbed by additive zero-mean Gaussian noise. 

The maximum likelihood estimation for the beacon coordinate- 
is found. The numerical technique used to solve the non-linear 
estimation problem is a modified Gauss-Newton method. This method 
consists of approximating the non-linear model with a linear 
expansion and iteratively solving the linear model using standard 
least squares estimation techniques. An approJiimate error co- 
variance matrix is found from which errors in the beacon coordi- 
nate can be estimated. Computer simulation has shown this 
technique to be stable for measurement errors within (and beyond) 
the range of errors encountered in practice. 

Many sources of error are present in the survey problem. 

The preceding sections detail the several steps to be taken in 
order to achieve an accuracy of ±1 meter in the beacon estimates. 
These can be generalized as follows: 

1. Sound velocity should be determined with a precision 
instrument, for example, a Conductivity-Temperature-Depth (CTD) 
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probe (see Gregg and Cox, 1971) . Errors in the measurement of 
sound velocity will usually bias the estimates but will not 
affect their variance. 

2. Corrections should be applied to the measurements to 
account for acoustic refraction. When these corrections are 
applied, the errors due to acoustic refraction can be consi- 
dered negligible. 

3. Travel time measurements should be made while the 
surveying platform is dead-in-the-water, or nearly so, to 
eliminate errors due to horizontal platform motion. This is 

a potential source of large errors in the pulse system, although 
it does not affect the accuracy of the Doppler system measurement. 

4. The time interval between selectecS Doppler survey points 
should be less than 1 hour since errors in the measurement of 
accumulated phase increase with time. Intervals of approximately 
1/2 hour each are recommended. 

5. Survey locations which will minimize the error magnifi- 
cation factors of the beacon coordinate estimates should be 
used. The choice of survey locations is critical to obtaining 
accurate estimates and can also affect the ability of the 
estimation technique to converge to a global rather than a 
local minimum. Several examples of good survey point selections 
are given in Section III. 3 above. 
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6. Pulse and Doppler measurements should be made simul- 
. taneously at no less than 10 survey locations. The use of 
'' more than 10 survey locations will, in general, reduce the 
I errors in the beacon estimates. For example, the use of 40 
' survey locations can reduce these errors by a factor of 1/2. 

In most cases the capacity of the processing computer will limit 
the maximum number of survey locations which can be used 
practically. 
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APPENDIX A 



DESCRIPTION OF PROGRAM ARCQR 

Pulse navigation systems calculate the distance between 
two points from the measured time of travel of an acoustic 
pulse between those two points. This distance can be approxi- 
mated by multiplying the travel time by an average sound 
velocity. This is an approximation since the sound travels 
in a refracted path rather than a straight line. The purpose 
of Program ARCOR is to improve the approximation by accounting 
for the refraction of sound in water. For a given sound 
velocity profile and the depths of a sound source and receiver, 
this program uses standard ray-tracing techniques (Officer, 
1958) to find the time of travel, TT, of an acoustic pulse 
between the two points is the travel time TT multiplied by the 
arithmetic mean sound velocity c: 

R = cTT . 

The actual distance, SR is given by the expression 

SR = (X^ + 

where X is the horizontal distance traveled and Z is the verti- 
cal distance between source and receiver. The ratio between 
the actual distance and the approximate distance, = SR^/Rj^, 
is calculated for several travel times. The set of approximate 
distances and their ratios (Rj^,Kj^), i = 1,N, are fit to a 
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Figure 7. Notation used for Program ARCOR. 
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second degree polynomial using standard least squares re- 
gression : 

SR 2 

K — — aQ + aj^R + ^2^ 

R 

where the a^ are the acoustic refraction coefficients. The 

actual distance SR can therefore be calculated when the 

approximate distance R = ^TT is known: 

2 

SR = R(3q + a^^R + a2R ) 

The set of approximate distances and their associated 
ratios are generated by calculating the necessary parameters 
for various angles 0 q at which sound leaves the source (see 
Fig. 7). For each angle 0 q there is a constant given by 
Snell ' s law 

sin0Q sin'0i 

= p (constant] = 

^0 ^i 

where c^ = sound velocity at the source and c. = sound velo- 
0 1 

city at interface i. This constant and the sound velocity 
gradients 




determine the travel time and the distance an acoustic signal 
must travel to reach the receiver depth, assuming constant 
gradient in each layer. The angle at which sound leaves a 
layer is 
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= sin ^ (pc) . 

The travel time for each layer is 

1 tan 0./2 

t = In ( ) . 

g tan 

The horizontal distance traveled is 

1 

X. = (cosG. , - cosG) (Officer, 1958). 

^ pg i 

For each starting angle 9^ the total travel time and horizontal 
distance traveled are found by summing the above quantities 
for each layer. We then have 

Total horizontal distance: X = ^ 

Vertical distance: Z = source depth-receiver depth 

Total travel time: TT = Xi t 

2 2 1/2 

Distance between source and receiver: SR = (^X + Z ) 

Approximate distance: R = "cTT 

Ratio: K = SR/R , 

The program computes the first starting angle by defining an 
angle, oc. = 90 ° and setting the departure angles equal to 
90° - c< . Each successive is reduced by a factor, for 
example .95, and a new departure angle is computed. The values 
of K and R are fitted with a second degree polynomial to pro- 
duce the coefficients, a^. A typical output is shown on page 

I 

73. A listing of the program begins on page 74. The pro- 
gram was originated by W.M. Max'quet of the Woods Hole Oceanographic 
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Institution and includes calculations of coefficients for the 
case of a receiver at the same depth of the source. This 
particular computation is not used in the survey problem and 
is not discussed here. The original version of this program 
was called Program SETUP and written by Mary Hunt and Roger 
Goldsmith of the Woods Hole Oceanographic Institution. 
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7.1564 


6 6.0 


i_*Zi0j'.ij.. 


._1 ■'‘8?.^l*7'b 


i 1:519 9. 't ■' 


111 Si-* :u 


7 . 3896 


6 7.54 


1*00013 


ll I 7 B . 54 



DtEP TO SHALLOW ( JO 5*d0 ) 

A SK ( 1 )_ r . 9^ 93 f 1 6 El 0 ^ 

ASK( P.i. : 5‘^0?^-90E-;-9 

A SK ( 3 ) = ,_2P1_ 3 Z 3.7E^ Ij 



Output of Program ARCOR. X=Horizontal distance between source and 
receiver;- SR' =' Slant range ; — rr="Travel time; - THSTA"= Angle of 
departure from source; K = SR/R; R = 2TT = approximate slant range. 
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FTN.L 



Oak’S 
OQk' A 



PROGRA M A RC OP 

DOUoi E PRECISION VERSION FOR [IP 

HP DO S VERS 1 0 M P,Y MARY KIMT & JAMFS [) U R [ l A H 



AUG 7't 



0 0‘’5 
BO” 6 



00k’ 7 
0 0" 9 

0 0] n 




FINDS RAY [-ENDliG COEFFICIENTS 

.DOUBLE. PREC ISI ON. PI t X I t.XF •.XD >J<MAX ».XC - Sk .PK •ThO.TH 
DOUBL E PRECISION THO^ * TH2 , T T » A RG > Al PH , AL PliT.RFAC 
DOUBl E...PREC 1 S.l OH XS»TS»CS^AL 2<RT 'ASt.tXT PACT »P, .sL 
dimension R'50' »DX(51') 

DIME N SIO N PTP C * 3) »QTSC I 

DIMENSIO’n Z.k 50 ) .C '5C ' 

DIME N S I 0 N _ G <_5 P > . __C BJ 5.B) 

dimension C0EF'3.A) 

dimension save '3) 

EOUlVALErjCi (CCEF; 1,]) >LTDC ) > 'COEF ( 

PI . 1 ‘‘ 1 5. ? G E 3a E il 



» 2 > >L TSC < . (C Oi F ( I > 3 ) ■ SA VE ) 



Nl 

NO 



SET UNIT 

: ] 

: T ' 



DEVICE COCES 



pause for sens:- SWITCFES 



^^PITE < fiO’ 

2 format !'• PROGRAM /RCO':« HP 

I" ssw I ON ‘ S'vU profile 

i” SSW 3 ON 

]•’ SSV; 7 ON 



yERSIOr;"''' 

LIST: L, or; 



SET 
UT F 



ifr.' 



. I T r t , 

r'-'iCF’ / 

■ FEEP TO deep CAL U'ATlCr; S^ i f F : L ’' / 

~ RAY [:EN[.ING CAL‘U[ATlCr. S )■ ; Pf E. I ' / 



H ^'!i 1 


2 


2 


C 


3 


C 




c 


' " 0 01* 5 




0 6 


1 ^'-5^ 



il^,_SS'A .Iji. or.’ - S. VE L PRO F lit 
v,P I Tt'k no'’ 1 1.^3) 
forma T(/" 'LL It. PUTS 
RE AD (tn IN!'’ ) IANS 
FCRMAT(-l) _ _ _ 

WRITE'NO'iai(t) 

_ E llM /: -T.< Ll ST DEV I CE” 
READ(r.I»^) l.P 



DO 2 I : 1 > 3 
COEF ( I » I ) : 1 .J: 

‘COEFLE’.lV ; 0.1* 

COEF ( 3 t I ) r 0 » k- 



ANC CC fl rs .C-UIP tl; .- 

PE free FIELL ’ /" r ‘ ADV ' ) 



At 



/ ) 






CONTII.UE 

INPUT SOUND velocity PROFILE 



l l“’ 7 
FF'* 8 
00 '*9 



WRITE 'NO'lPk'5) 
f .ORMAT (’• S._. Y.Fl 

READtr-I.*) I SVC 
wRITE'NO’Il'-'b) 

FORMATt" ir.'PuT PRCflLE"^" 



INPUT DEVICE? -;■) 



LR PTIK f ) 



. vr L rr /^. 
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1" ENTER NEGATIVE NUMtiER VvHEN D0NE!"> 

li 0? 1 I : P 

00 ^ 2 '" 5 cbNTItiUE 

0 0^3 I = ! + 2 

OfiV‘1 READ ( ISVC » >* > 2(I)’C(I) 

10^_5.. IF |2 ( I ) LJ:218i.fi ^ 

0 0^6 8 CONTI t;UE . 

0 8*^ 7 C : ^ 

0Cr‘8 C no’ MORE THAN 50 ALLOWED 

0 0^9 C 

0Of-r: IF (I - 50) ‘-’»3'-'-.3P 

12_C0.NTIt;uE_ 

00 t ’2 ’ 1 = 1^1 

£Bt3__C 

0 0‘''4 C 

0095 C IN ITIA LT 2E VARI/PL ES 

0096 C SET X INITIAL’ TERMINAL AND 1 L;C REMFf,T aL VALUE 

0 09- 7 C 

00< 8 30 CONTir.UE 

0 0i^9 NS V : I _ 

00/ ,’” ' IF ‘ ISSW ( 1‘) ) 31 ’ 32 ‘32 

00/1 31 WRI TE ' LP’30 00) 

00/2 3000 FORm/'t’l/:' I 2 C’/) 

S 0/ 3 ___ wR I TE ‘ Lf^’_3EPJl ( I’2II)>Ci:)»I =l » r;SV ’ 

’ 0 ’ 0 /it ’ 3-"’0T format i r- > 2 ri ’0 *‘L ) 

__ 32 IF'LSSWC,’)) 2 01-T ,.3.'t , 

00/6 3<4 COt.TlEUf.' 

0 0/ 7 XI ; LL’-LDP 

00 / 8 XF -‘ 200 0 • OS' 

^0^9 _ xr = 2090 . DP 

009 MORD : 2 

0 0 t- 1 _ _ MORDl =_MORC_ + _1 

000 2 ”5G ’ continue’ 

F'00 3 IT ER - “1 

« 08 'l C IL.PUT typical DFlP DEPTH 

0095 WRIT^'NO_’^^ 0 •‘' 0 ) _ _ 

e'006 ’ r’000 FORMAtV" beacon CEP’fh^ -"’) 

0 00_7 RE AD (N I.’) 2 0 

0008 C ir.PUT StiALLCW RECEIVER DEPTH 

0 00 9 wRITE'tJO’ 4010) 

009.-' /) 0 iB format (" depth of RECEIVlR 1? ) 

£091 READ (Nl 2R 

0 0 V 2 " " C 

00 M 3 c +’+ STARJ deep TO DLIE P__PRqCE5D5 1 NG 

’009 'I ’ C ’ • ” ” 

00 M 5 Q find AVERAC.E ViLOCITY AT TYPICAL PFEP DEPT [1_ _ 

0 0 ” 6 7 0 COtJT’lLUE 

0P’'7 IF'rSSwL’)) iri'0»75,‘‘D 

0C'”8 75 CONTIiTJE 



-75- 



■ on 




: 

I 

: 



f 



I 



I 



page HB03 ArCOK xnos - FORTRAN IV COPf’UFP 



00'^g CPAR = Yp(Z0.ZiC»NSV' 

0 1 » i i c P I N D _j hq:_ d r ep depth starting l me R 

0101 bo 110 II : 1*NSV 

01 02 IF ( Z0 - Z ( 1 1 ) ' 1 1 5 » 1 15 vU0 

0103 110 COfjTIf.UE 

0104 i.L.^tLsy , 

0105 ?15 CONTit.UE 

011V6 C 

0107 C find gradient of sound velocity at deep depth 

010 8 Gdl) = 'C(II) - Cdl-lM/'Zdl) - Ztll -lM 

0109 c 

01d] C__ IF__GRAniENT_JS__K'EGATI VE*_DEEP_WATER_ANAI YSIS_iS_ SKIPPED 

Pill IF (Gdid 125'130 .i35 

_^ll2 ^1.2 5_ WRITE. (MO <9.1.0 2 > 

PiTb 9102 format(5'h (jegatIve gradient ~ DEEP TO deep processing' RYPASSEC) 

0 11 G O TO 10 -'P 

PII5 130 WRITE (NO<9101' 

EJ;6 9}21 FCR.HAT(92p GRACIEMT_: P FOR DEEP TO DEEP PR0CESDSING7 . 

p'li? 20H K : l.C 'for ALl X >' 

Bl£8 GO TO 10 0.0 _ _ 

011 9" C 

0 1^ C ( A LC U 1. A T E__E 0 L CO P__ P PE TERS 

'PIE-I C 

jnE2 _ 1_35 CONTir.UE 

0iE3 lE : ,‘XF - XI )'/xO +~’0.0. ID" " 

jdl'L- c SET CO I COP TO Vary x for deep to deep 

fdVs WRITE (I PtQins’ ■ 

9’. -iS FORMAT (//// ‘- X 1 P X 1 '' K " > 1 H X K 1 ' 

'"Db''4S;’“r ="l'dc 

01-"8 FI : I _ 

rJlb'y ■ ■ xc' : XI E'fI'XC ’ 

•d-V IF (XC) 17Sd'-'d75 

0di>i C SET THETAC FOR X ■= 0 

iiiP2 150 C OM Tit ^ ; 

01:’ 3 DK = -1 . PDP 

01:’<t_ Rd) = 0.P _____ ___ 

Ini’S ’ " ' GO TO 35P 



PJA-G ^ CALfUlATE THETA0 

01? 7 ” 175 CONTIdJE 

01r’H TPi- : 2 . 0 < C CAR A ( XC ^ G d I > ) 

01:’9 TI-10 r DATANdHP) 

0 i‘d c calculate angle TH 

Pi:* 1 - 20 0 CONTIrjUE 

01“' 2 TH r PI - TH0 

" 0 IR 3 C FIND half angles 

01“ 4 ■ TH02 ~ TH0/2.0r (! 



5 


TMi 


2 : TH/2.0D0 




C 




calculate time t 






TT 


: .'DSlfdTIiP’/rCOS.dl:*') ) / .'DSIN(THc 


<•) /( COS ( Th 2> ) 
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01i'8 

01^9 


c 


TT r DLOGCTT) /G( II ) 

CALCULATE APPARENT DISTANCE TRAVELLED 






RACT = CDARm 
R(T) = RACT 


01^3 


c 


calculate RATIO K 
DK = XC/RACT 


01^5 


•550 


CONTIt:UE 

DX > I ) = DK - 1 .0D0 


Hl?6 

f)l^'7 


•'0 0 


WRITE (LP»9i06l XC»CK»rx(I' 
CONTI DUE 


in? 8 
13ir)9 


c 


FIT polynomial with X AND K TO GET COEFFICIENTS 
call HCRFT(R»DX>r,TDC»lE -MORD) 


filo ,, 
ijjtU 


c 


OUTPUT CURVE FIT COEFFICIENTS 

DTDC (3 ) - DTOC <1)41.0 _ 


01(^2 
0 3^ 3 


c 


WRITE (LP»9l03> Zk- , ( I ,DTOC ' I ) ' I = I'MORDl' 


01f'9 
r ] o 5 


c 

c 


♦♦ DEEP TO shallow CALCULATIONS 


6 

fM^’7 


1 01; 0 


CONTli UE 
CS : 0 . or 3 


V ] '> 8 

y 1 1*. 9 


' C 
C 


FIND STARlir.G DEPTH AND LAYER 

Flr.’D layer GfiADIEf.T AND AVEFA&f SOUL.D VELOCITY 


*-* 1/ • 
mil 
*‘1^2 
Wli 3 

>n m 

m/5 


1 n 0 
1 i 5:^ 


DO 1120 j; 1 ,NSV 

IF (Z0 - 7(J)) 11E0 ,1 1 8i' I 11 00 _ 

CCr.TIL.UE 
J - HSV ^ 1 

COL.TIi UL 
K - J 


0 1/6 
01/7 




ZJl : ZLJ-I' 

Cr-AR - YF (Z0iZ j>C .NSV > 


6] /s' 

0 1 / 9 




ZJS r ’Z(J) 

CJS r C ( J) 


iU< , 
aibi 




JS z J 

Z(J)=Z0 „ 


0 1 1> 2 
1 1 3 




C(J) : CLiAR 

G(J) = (CDAR - C (J-1 ' )/ :Z( J) - ZJl) 


0 1 ^ 4 




CFr'J) z' 'CPAR •• C'J-1))C2.0 

GO TO 121'0 


01^6 

01^7 


c 

1 180 


FIND layer gradient AND AVERAGE VELOCITY 
JS z J 


016 8 
f>169 


1200 
120 5 


ZJl z Z<J-1> 

CONTli. UE 


m''' .: 
01‘^1 
0 1 2 
file's 


120 a 
1 2] 0 


K :: J 

IF (2'J) - ZR) 1300il30i’4lZ0B . 

IF (ZJl - 'ZR) ' 3 210 . li-’30 ’ 123 0 

CONTiruE - 


0 J ' ' 4 
Pl‘^5 
01^' 6 




CD'AR : YP (ZR *Z iC «LjSV ) 

G(J) = (C(J“ - CHAH)C(Z'J) - ZR) 

CP ' J) z 'CP AR * C ' J) > /2* 0 . 
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K : 1 

ZJl = ZR 

mv9 GO TO 125 15 

02''^ O' GU) = (C(J> - C(J-JM/<Z(J) - ZJl) 

02'' 1 cn' J) = ‘C(J) ^ CiJ-3)) /2.t' 

e[2J12 1 ^50_CONTIt:UE 

020 3 ’ C5 = CS CD(J'4(Z(J) - ZJl) 

0204 _ 

02^5 ' " j' :' K ' "" " 

0206 I F *J “1) 130 0 >1275«12'00 

020? 1075 ZJl : -O*.0.'t: 1 

^±8 G 0_ T 0 12 05 

02^9 1-500 C5’;'CS/.'Z0 -ZR)' 

02.1 j3 WRITE'LP >910_9) CS 

’0'2ii 'i’-5'2'e contii-ue 

02j 2 IT ER - ITE P 1 

0213 C SET STAR-ilrjG Af.'GlE “ 0 DEG 

0214 At PH ■ 0"Or. !' 

0?i 5"” ' Al'PHT : PI/2.CO0 

02)6 _C _ SET ARRAY ItJEEX 

021 7 1 = 1 

0?I 8 BI.A.Q J . - ^ ‘ 5 c ^ 

02? 9 ' WRl’fE~TJO’9201) ' ' 

g2r-l PORVATI" ENTER CFCP^PENT FACTOR BETvs-'Frj 0 ANT 1 

?-2c'l READ (fa 1 ’ ) RFAC 

PZ'r'2 . __ 'aRiTE ‘ f;0’92R2) _ 

’02E3 9P'.'2 FOPHATC EfT.R M A X HORIZONIAl RAt.cr -".) 

i:^4 . XHa:. . = IJNE.V.L- 

R2-E5 read ua»*) XPAX 

R2t 6 WRITE' LR’ 91 7 ) 

r2i'7 C 'LOir. RAY trace 

0 2E8 1“'? COfTI, UE _ 

02E9 C IiaTlAlizE SUES 

0?^ XS_= i’.0C? 

02:5 1 TS - l‘.0C0 

P-2.5 2 S - p.l7r“ 

02-‘'3 PrDSif. t ALPH f /C ' JS > 

0254 Tl-TTA =._ALPH + 1 E0 • 0D0/PI 

02;55 J : JS 

0 2~-6 K = J 

02^7 1‘-20'ZJ1 : Z(J-1' 

0258 l‘>25 COtJTH.UE 

0 25 9 CJl :CtJ"l' 

02“.. IF ‘Z(J) - ZR) 15rn5.lNfi;-,,’1423 _ 

02'Jl 1 420 IF 'ZJI - ZR ' 1 430 ’ 14"0. I'HH* 

Jl24 2 1 4 3 0 COfJlI.L UE 

024 3 K = 1 

0?"4 ZJl = ZR 

s'24 5 CJl = U-AR 
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»21*6 
UpM 7 


l'»40 

C 


COrJTIfiUE 

find angle ray leaves layer 




2:.» 8 
9 




ASn = Cjl+P 

TF (ASM - 1*flD2) ] 490, 1 'Itspn 




(i 2^ 

8 2^7 1 


145C 

C 


continue 

TH : AS IN (ASN' 




J5 2^2 
2*' 3 




ARG : ASN/DSORT ( 1 -GDi'- - ASN+ASN) 
TH r DATAN(.'.RG> 




B2? 4 
112^-5 


i‘‘6n 


IF'G(J)) 1462*1480,1462 
CONTI E’ur 




i5 2^>6 
02^7 


c 


find TU'E 
TH2 r TH/2.CD0 




02^8 

{'<279 




AL2 : ALPH/2.0r;0 
IF'1“1) 14 72*1400*1 '(70 




1>2^' 

.^2<'l 


1 472 
C 


C0t;TH.UE 

TT : D7Ar;(TH2)/DTAN(/L2’ 




H2'-2 
i’ 2*' 3 
i)2^' 4 
02t'5 
6 

i'2^’ 7 




TT : 'D5IN(Th2'/LCOS,'TH2) ) / 'DSIN( 

T T = ' -n 06 fjT > 1/G(J 7 

RT t ' -1.2D4) /*p4C-(J' ) 

XT r 'DCOS(ALPh) -[<COS(T|<) ) 'RT 
AL = < TH'Al Ph) 'RT 

GO T' 14‘'< 


AL2' /LCOS ' AL2) ) 


0 2^-8 
P 2 f'' 9 


c 

^ U p a 


calculate L-AYF.R values FOR STRA 
ccr:TivuE 


IGM UP OR G r CASE 9 


0 2/ 1 
!^2/2 " 
r- 2 / 3 


c 


XT : '2 ( J)- Z *DT At; (■ L PH ' 

XT = 'Z(J)- ZU! ) * iDSlN(ALPh) /ECOStA 

Al V XT " ' 

TT = lSORT '(Z'J)- ZjDMZ'J)- Zul) 


IPH') 

' XT ( - ; ) Ac ! ' 0 ) 


;■» 2i 4 
0 2/ 5 


c 

1<‘9£ 


CONTI;. UE 




' .';276"” 
‘'2/ 7 


C 


SUM X EIS' TAfJCE 

XS : XS ■* DABS' XT) 




r<2/ 0 

i)2/ 9 


C 


SUM TRAVEL Tl"F 
TS : TS ' DABS'TT' 




0 2V' ■ 
2^ 1 


C 


SUM ARC LEi;GTH 

S : S ^ AL 




0 2^ 2 
0 2^^ 3 




RESET Entry PAY for next layer 

ALPH - TH 




8 2V 5 


T 


re’set index And check 

K r K-l 




{i29 6 
f! 2«< 7 




U : K 

IF (J - 3) 3 501' * 1496 ’ 1 4 20 




fi'2^'0 ’ 
1)2!; 9 


’ 1495 


Z.Jl r -P- F 
GO TO 1425 




f/2"'2 


c . 

1 ‘='0 0 


STORE totals 
CONTI t.UF 




0 2‘^2 
02'*3 
i-2'-'4 




R A C T - X S ' )• S ' ( Z - Z R ) t • Z P - Z R ) 

PACT : rSCPTtMACT) 

R ( 1 ) : C S ' T S 
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H2‘^5 

82*^6 



SK = RACT/RtI) 

_WRITE ‘ LP’O!'"?) XS’RACT »TSiTHETA .SK ,R ( I < 



02V? 9 IB 8 FORMAT (Pl 0 . .n C .2 , Fl fi.'t ,f 1 0. 2 ,F 1 fi .S ,F 1 0.2 ) 
H2^8 DX < I ) = SK ~ 1 . 0PPi 



0 2^9 
03*^ 0 



j 

0 3^^ 2 


0 32 3 
0 32 4 


1E5C 1 


0 3i^5 


c 


0 3^^ 6 




03^7 




03.2..8_ 


^ _ .. 


0329 


It 5 0 1 


03 i 


1 


03d 1 


c 


^5 31 2 


j 


i>3i 3 




0 31 






y/vi' < 


6 


. 1 


M 3 i 7 


( 


^3_i_8 


25 


i:' 31 9 


1 


T: 3^: 


i^O2 0 1 


• 3r j 


i ' 


: f ? 


1 ' 


bZt 3 


...... 


0 3- ^ 




5-3£'3 


195 0 I 


t 3l 6 


1 


0 3‘^"7 “ 


1 


i* 3*'' 8 


1 


0 3^ 9 


2t*&0 


(•'. ~t Si 




03 -1 


”2*'10 ' 


H 3 :'^ 2 




■'03<^3 


' 21*12 ' 


0 3:=>^ 




0 3-^5 




0 3^^ 6 


2»^/l 


0 3-^7 


1 


0 3:^ 8 


1 


03^9 


”5000 


0 3'^ 




0 3** 1 




0 3“ 2 




0 3‘^ 3 


5 ‘*4? 



REINITIALIZE AMGlC AND INDEX 

=_ J „t_i 

ALPhT : RFAC4ALPHT 

Al PH : PI /2-0D0 - AlPi-lT 

IF.'XS - XMAX) iSSU 

COfiTIt.UE 



CHE^K that INDEX DOES NOT BECOME TOO GREAT 

IF (I - ±Q I , 

I - I 



I : 

F IT P.Ol YNOHI/ L VnT!i_X_ANp_K_£OEF.Fl CI.FNTS_ 
COtJTlfUE 

CALL J!CEi'.T.,lRj.rX„.iJSC.'J iL-jDPX.) 



OUTPUl CURVE FIT COEFFICIENTS 

DTsc(i) - rTsc'n -1 _ ■_ 

WRITE (Ip.qlfiij) Zi' .ZR , ( I I ,DTSC ( 1 I ) . II .MORr I ) 

IF < ITER) 58’I,7''f;,2 IH’ ... _ 

coririi.uE 

.Dp 2 S I 

SAVE ( I ) : CTSC ' I ) 

COt.T I' UE 

WRITE 'NO. '122 n) 

FORMAT (/'■ depth OF RECElVI. R 2^"/ 

" 'EDI. R r.EOAVIVi NUHCDR 1( .-NOIhEF set OF' 

. c Cl f,F 5..I „..f!p.T . Ft-.siRD F ' 

readini.') zr 

IF ‘ ZR ' I 9f)>' . 1 DO . 1 : i.' 

COt.TIiUi 

DTSC n ) = 1 • 2 

CiTSC (> )' ■ 'd *2 

pi SC(T ) = ('» £ 

CONTI'.UE' 

IF t i SSw (U|) » 2i'10. 7: Pfi ..7_..pD 

continue 

If ‘ZP - ZtJS)) 2Pi 4 .2Pl^.2Fl't 

COfjT If UE 

Z(JS) r ZJ5 



C ( JS) = CJS 
CONTir:UE 
NPT r 'I 

WRlTE'NPT.5i'00) NSV 

format (IE) 

DO_'l5_ I =1..J1SV 

C ( I ) : C ' I ) / 1 B 0 . 

WR] TF ' NPT .5i'4F: ' 7 ' I ) ’C ( 1 ) 

FORMAT (Dr 1 C.'( ' 
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0 3^ 4 
C3*‘ 5 


'15 CONTIflUE 

5060 .format ( 3E20 • IjlL 






0 3'* 6 
0 3'* 7 


WRITE‘nPT»5060I 
WRITE ' NPT .5060 > 


(SAVE(j' .J=1.3) 

!nTSC( J' .J:i.A) 





03^8 write ‘NPT I ( DTDC ( J ' i J-1 . ’ ) 

0 3‘< 9 7 0n0_COtlTI!iUE_. 

0 3‘Jf; STOP 



0 3" 1 9,LB;3 . FORMA Tjy//2'^H RF,.EP_ T 0_DEEP_ _Df.PTH 

0 3^2 ‘5M AOK ( » I2,3tH) 

03^3 Q l.R'j FORMAT < /// ?IH D EE P To Sl- ALLOW ■ iF10-2 «3H TO'EH'-RiFh >/ 

f)3P0 ^ (SH ASK ( ’ I2> 3H) =>E1'I.7M 

B_32_5 99 E 6 fO R M A T ( 3 E 1 2 . 5 )_ . 

0 3" 6 9^E7 FORMAT ( // /5X , 1 HX t iSR » 9X > I HT ^ShTHETA ^ 5X - Ih-K . PX > ' HR ' 

.0317 9309 FORMATT /" arithmetic PEAN..SOU(lD, VELOC ITT : " , ri?»A).., 

C3"'8 EMD 
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0 0*^1 


c 




curve-fitting SUE’ROUTIME F*K.KEYTr UN‘H.0. 


I • 


) 


12TH MAY l98‘l 


P 0 V ’2 






SUBROUTINE IICRFT ( X » Y ’ C . N » MCURV ) 








0003 


c 




VERSION FOR ItEWLl TT " PACKARD 








D 0P /( 


c 




THIS IS A Cut-down version cf suiiroutinf crvft 








0 0^-15 


c 




DONE TO SAVE ROOM 








00»^6 


c 




DONE BY M. HUNT MARCH 19 / ^ 








001 ’ 7 






dimension .X I 1 ) » Y ( 1 ) 1 C 1 1 ' 








0 0 I-' 8 






DI mens I on a ' 5 1 (• ) . AP ( s » 6 > 








001 ’ 9 


c 




M IS COEFF INCCTR’ 1=1 IS SPCL ROUTE INDCTR FOR 


Y = 


cdi only 
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APPENDIX B 



DETECTION PROBABILITIES OF WOODS HOLE OCEANOGRAPHIC INSTITUTION 
CYCLE COUNTER 



The detection probabilities of the Woods Hole Oceanographic 
Institution cycle counter can be calculated by considering the 
problem of detecting a sine wave in the presence of Gaussian 
noise. This treatment is sufficient because the estimate of 
phase is not based on the optimum estimator 

-1 

0 = tan 



Si (t) 

where Sj^(t) and S 2 (t) are the quadrature components of the 
input to the phase detector; the estimate of phase is based 
on the signs of s^(t) and S 2 (t) . The detection probabilities 
of the sine component will be developed for illustration. 

The corresponding probabilities of the cosine component can 
be developed in a similar manner with identical results. 
Consider a sine wave in the absence of noise 
s (t) = A (t) sin 0( t) . ■ 

The amplitude A(t) varies slowly compared with the phase 0 {t) 
and can be considered a constant. The phase can be assumed 
uniformly distributed 

p(0) = 1/2 -rr ,\ 0 \ < TT 
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V) 



The probability density of the signal p(s) is therefore given 
by the expression; 




The phase detector output in the presence of noise will be 



denoted by z(t) = s(t) + n(t) where the noise is normally 



distributed: 




-n^/2 O' 
n 



2 



(O = r.m.s. noise voltage. 



The probability density p(z) may be found by the convolu- 
tion of p(s) and p(n) assuming statistical independence of the 
two processes. This leads to the following integral for p(z): 



p(z) = 



1 r -r?/2 0'^ 1 

n 



n 



J 



’fTA 






^ (z-n)^ 



This integral has several representations (se'e Ol'shevskii, 
1967) but is approximately a constant near z = 0 for ^ 1* 



Figure 8a illustrates the shape of this probability density 

2 

as a function of the parameter = "2q^ * 



The probability of error is a weighted sum of the proba- 
bilities of false alarm P(FA) and false dismissal P(FD). The 
cycle counter decides that the signal is positive when the 
voltage exceeds some threshold > 0.. Similarly the decision 
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that the signal is negative occurs when the voltage is less 
than some negative threshold -Vrp < 0. A "false alarm" occurs 
when the signal is below the threshold and the signal plus 
noise exceeds the thresholds: 

|2(t)| > given |s(t)| < . 

The probability of this event, P(FA), is approximately the same 
as the probability that the noise will exceed the thresholds 
in the absence of signal. This is represented by the shaded 
area in Figure 8b and is given by the expression 

r r ^ 

P(FA) Cl 2 / p(n)dn = 2 exp {-x^/2) dx 



V. 






A false dismissal occurs when the signal is above the 
thresholds but the signal plus noise is below the thresholds: 
jz(t)| < Vrp given s(t) > V.p. The probability of this event 
P(FD), is approximately the same as the probability that the 
signal plus noise is below the thresholds for any signal. This 
is represented by the shaded area in Figure 8c and given 
by the expression 



V„ 



P(FD) 



p (z)az . 



-v_ 



In the region of interest p(z) can be approximated by a constant 

1 



P(z) 



TTA 
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Figure 8- Probability Densities of Signal Plus Noise 
p(z) and noise p(n). 
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The probability of false dismissal can then be found 

2 2 (V^/cTn) 

P(FD)2: ''T/ = ^ . 

-tr (A/ O' ) 

' n' 

It can be seen that an increase in the threshold voltage- 
to-noise ratio Vip/ O' ^ will increase the probability of false 
dismissal and decrease the probability of false alarm for a 

given signal-to-noise ratio. Since the losses assigned to 

I 

these errors are approximately equal, the optimum threshold 
setting for a given signal-to-noise ratio will satisfy the 
following equation: 



P (FA) 
P(FD) 

or equivalently ^ 



= 1 



1/(2-Tt) e • dx = 



V, 



T/C" n 



V, 






■n (A/q;) 



( 10 ) 



The expression on the left hand side of equation (10) is a 

tabulated function of V,p/ and can be solved by trial-and- 

error for various values of the signal power-to-noise ratio 
2 

S/N = A^/2 Cf^ (see Table 7) . 

The probability of error P^ in determining the sign of the 
signal is given by the following expression for the optimum 
threshold setting; 

Pg = 1/2 [P(FD) + P(FA)] = P(FD). (11) 
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TABLE 7 



10 log S/n 


''vcr„ 


P 

e 




12 


1.40 


.158 


.366 


18 


1.67 


.095 


.208 


24 


1.92 


.055 


.116 


30 


2.16 


.031 


.064 


— * 









An error in the cycle count, i.e. the number of quadrants 

of phase change between two successive cycles can be caused by 

an error in either or both of the quadrature components. When 

an error (false alarm or false dismissal) occurs in only one 

component the cycle count will be in error by ±1 quadrant. 

When an error occurs in both components the cycle count will be 

in error by ±2 quadrants. The following probability of cycle 

count error can be calculated from the probability of error 

previously determined by equation (11)1 

P(e=0) = (1-Pg)^ 

P(e=±l) = 2P^(1-Pg) 

2 

P(e=±2) = P^ 

The variance in the n^^ sample can be calculated as follows 
^ = P(e=0) (0)^ + P(e=±l) (1)^ + P(e=±2) (2)^ 

= 2Pg(l-P^) + 4Pe^ = 2Pe(l+Pg). 

The variance of error in the measurement of accumulated phase 
change 0rj> after T seconds is 

= E ^ = N " N = Tf^ 

•^T n=l n 

where fg is the sampling frequency. For a sampling frequency 
of 10 Hz and an interval of 1 hour =’3600 seconds the standard 
error in the measurement of 0^^ is z 190 , 
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Table 7 lists the values of P , ^nt ^ and the optimum 

e 

threshold-to-noise ratios for various signal-to-noise ratios 
typically encountered in practice. The variance is calculated 
assuming that the bias is zero or negligible. This assumption 
is valid when the velocity relative to a beacon is small. 
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APPENDIX C 



SAMPLE CALCULATION OF SIGNAL-TO-NOISE RATIO 

In general the signal-to-noise level will decrease as the 
range between a sound source (e.g. CW beacon) and a receiver 
(e.g. ship transducer) increases. A sample calculation is 
made to demonstrate the realtionship between signal-to-noise 
ratio at the source (S/N)g^the range between source and receiver 
R, and signal-to-noise ratio at the receiver, (S/N)^. The 
signal level of the CW beacons in the present configuration of 
the pulse-Doppler system is 166 dB (re: l^Pa). A typical 
ambient sea noise level in sea state 3 is 44 dB/Hz (re: 1 Pa) 
or 54 dB/lO Hz (see Urick, 1967). The signal-to-noise level in 
a 10 Hz band is therefore 112 dB. Transmission losses are 
assumed due to spherical spreading and attenuation. The loss 
due to spherical spreading is 20 log R where R is the range in 
meters. The attenuation of a 13 kHz signal is approximately 
1 dB/km. The signal-to-noise level at the receiver (S/N)^. is 
therefore: 10 log(S/N)j. = 112 dB - 20' log R - R x 10 ^ . For 
a horizontal range of 10 km and water depth of 5 km the range 
between a bottom-moored source and a surface receiver is 
approximately 
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R = (5^ + X 10^ = 11.2 X 10^ m 

The (S/n)^ in decibels is therefore 

10 log(S/N)j- = 112 - 81.0 - 11.2 = 19.8 dB. 
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